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Stability is a desirable property of complex ecosystems. If a community of interacting species 
is at a stable equilibrium point then it is able to withstand small perturbations to component 
species’ abundances without suffering adverse effects. In ecology, the Jacobian matrix evaluated at 
an equilibrium point is known as the community matrix, which describes the population dynamics 
of interacting species. A system’s asymptotic short- and long-term behaviour can be determined 
from eigenvalues derived from the community matrix. Here we use results from the theory of 
pseudospectra to describe intermediate, transient dynamics. We first recover the established result 
that the transition from stable to unstable dynamics includes a region of ‘transient instability’, where 
the effect of a small perturbation to species’ abundances—to the population vector—is amplified 
before ultimately decaying. Then we show that the shift from stability to transient instability can be 
affected by uncertainty in, or small changes to, entries in the community matrix, and determine lower 
and upper bounds to the maximum amplitude of perturbations to the population vector. Of five 
different types of community matrix, we find that amplification is least severe when predator-prey 
interactions dominate. This analysis is relevant to other systems whose dynamics can be expressed 
in terms of the Jacobian matrix. Our results will lead to improved understanding of how multiple 
perturbations to a complex system may irrecoverably break stability. 


I. INTRODUCTION 

From the perspective of local stability analysis, if an 
ecosystem is close to a stable equilibrium point then the 
effect of a small perturbation, such as the loss of indi¬ 
viduals from a population, will eventually decay and the 
system will return to its original equilibrium point mm- 
But if the ecosystem is at an unstable equilibrium point 
then the perturbation will lead to the system settling 
at a new equilibrium point, possibly with fewer individ¬ 
uals or even species mm- In theory, ecosystems with 
large numbers of species and interactions are more difh- 
cult to stabilise [5]. However, many ecosystems contain 
vast biodiversity mm- Reconciling this finding with lo¬ 
cal stability analysis has motivated ecologists for over 40 
years [S]. 

Recently, stability criteria were extended from ran¬ 
domly assembled communities to include those with more 
realistic compositions of mutualistic, competitive and 
predator-prey interactions [5] . These criteria showed that 
communities in which predator-prey interactions domi¬ 
nate are more likely to be stable. It was then shown, 
using empirical food webs, that the distribution and 
correlation of interaction strengths has a greater effect 
on stability than topology: how species interact with 
one another is more important than who they interact 

with nnun]. 

Stability is a long-term concept: it indicates whether 
a system will, at some point in the future, return to the 
same state as before a perturbation [T^- Reactivity, on 
the other hand, indicates how a system will respond im¬ 
mediately after a perturbation has been applied [IMZ]. 


A stable system can be non-reactive, meaning that a per¬ 
turbation to species’ abundances dies down immediately, 
or reactive, meaning that a perturbation is first amplified 
before eventually decaying (whether a particular pertur¬ 
bation is amplified in practice depends on which species 
are perturbed and by how much m)- Reactivity crite¬ 
ria for large ecosystems show that communities on the 
verge of instability exhibit reactive dynamics CHI: and 
identifying a system as reactive has been proposed as an 
early-warning signal for population collapse [T2H23]- 

The starting point for deriving criteria for both stabil¬ 
ity and reactivity is the community matrix |24j . A spec¬ 
tral decomposition of the community matrix provides in¬ 
formation on the asymptotic behaviour of the system for 
stability {t —)• oo) and reactivity (t —)■ 0). But so far, lit¬ 
tle information has been extracted from the community 
matrix regarding transient dynamics: how the system 
evolves after a perturbation and before it either returns 
to equilibrium or becomes unstable 

Reactive dynamics are not possible if the community 
matrix M is normal, i.e., MM^ = MlM, where is 
the adjoint of M [2S]. But if M is a non-normal ma¬ 
trix, as is usually the case in analyses of realistic ecosys¬ 
tems, then transient dynamics may substantially differ 
from the asymptotic behaviour suggested by the eigen¬ 
values of M. In addition, small changes to the entries 
of non-normal M can cause an otherwise stable matrix 
to become unstable [IS]. In such cases, the dynamics 
implied by non-normal matrices are better described by 
pseudospectra, which detail the neighbourhood of eigen¬ 
values in the complex plane for different average changes 
to the entries in M [21] . 
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Here we formalise the transition from stability to insta¬ 
bility in terms of pseudospectra. Using this approach, we 
consider the effect on dynamics of two kinds of perturba¬ 
tion: more commonly studied perturbations to the equi¬ 
librium abundance of species (to the population vector) 
and less commonly studied perturbations to the entries in 
M (which could be interpreted as uncertainty in, or small 
changes to, species’ interaction strengths [30]). We de¬ 
scribe critical values for community properties separating 
three regimes: stable and non-reactive dynamics, stable 
and reactive dynamics—‘transient instability’—and un¬ 
stable dynamics. We show that system dynamics at the 
boundary between non-reactive stability and transient in¬ 
stability can be affected by perturbations to entries of 
the community matrix. And, given a perturbation to 
the equilibrium abundance of species, we provide upper 
and lower bounds to the maximum amplification of such 
perturbations during transient instability. This allows us 
to sketch out the transient dynamics of complex ecosys¬ 
tems using only information from the community matrix. 
Finally, we compare the properties of community matri¬ 
ces representing ecological communities with five differ¬ 
ent types of interaction structure: random, mutualism, 
competition, mixture of mutualism and competition, and 
predator-prey. 


II. METHODS 


A. Local stability analysis 


Here we consider an ecological community of S species 
for which their population densities at time t are given 
by the vector Y(t), as in Tang & Allesina [H]. The 
dynamics of the population vector Y can be described 
by a system of coupled differential equations 

f = <» 

where f = [/i, /2 • • • , /s]^ is a vector of linear or nonlin¬ 
ear functions. An ecologically-relevant equilibrium point 
is a non-negative vector Y* such that 

f(Y*) = 0 (2) 


The community matrix M is defined as 


Mij — 


dfz 

dYj 


Y=Y* 


( 3 ) 


which is the Jacobian matrix evaluated at an equilibrium 
point |24| . It is well known that an equilibrium point is 
(locally and asymptotically) stable if any infinitesimally 
small deviation, AY(0), eventually decays to zero, i.e., 
limt_>oo ay (t) = 0 125 . In the vicinity of an equilib¬ 
rium point, the time evolution of a perturbation can be 
described by 


( 4 ) 


Therefore, the spectrum of the community matrix M is 
clearly relevant for determining local stability. If A(M) 
is the set of eigenvalues of M, then an equilibrium point 
is stable if all eigenvalues have negative real part, i.e., 
Re{X) <0VAe A(M) 011]. 


B. Generative models for community matrices 


We parameterise community matrices using four quan¬ 
tities: S, C, /r and a; where S, as above, is the num¬ 
ber of species, C is the connectance (the fraction of re¬ 
alised interactions among species), /r is the strength of in¬ 
traspecific interactions and cr is the standard deviation of 
the strength of interspecific interactions [1]. We assume 
that populations are self-regulating and so Mu = — /i, 
where fi > 0. Non-normal community matrices with dif¬ 
ferent types of interaction—representing different types 
of ecological community—are generated by sampling off- 
diagonal entries {Mij , interspecific interactions) from dif¬ 
ferent bivariate distributions. Having specified a partic¬ 
ular distribution, stability criteria can be expressed in 
terms of S, C, fj, and a. Based on these criteria, it has 
been shown that predator-prey community matrices are 
the most stable, followed by random, competition, mix¬ 
ture and mutualism |1] . Generative models for these com¬ 
munity matrices are described below. 

Random. Each off-diagonal entry is sampled indepen¬ 
dently from a normal distribution ^(0, a) with probabil¬ 
ity C, and otherwise Mij = 0 with probability 1 — C. 

Mutualism. Each off-diagonal pair {Mij,Mji) is sam¬ 
pled from a half-normal distribution 1^(0, cr)| with prob¬ 
ability C, and both entries are zero otherwise. These 
community matrices have a (-I-, -I-) sign structure for off- 
diagonal pairs. 

Competition. Each off-diagonal pair {Mij, Mji) is sam¬ 
pled from a half-normal distribution — |A/'(0, (t)| with 
probability C, and both entries are zero otherwise. These 
community matrices have a (—, —) sign structure for off- 
diagonal pairs. 

Mixture of mutualism and competition. Each off- 
diagonal pair {Mij, Mji) is sampled from a half-normal 
distribution 1^(0, (t)| with probability C/2 or — |A/’(0, a)\ 
with probability C/2, and both entries are zero otherwise. 
These community matrices have a (-l-,-l-) or (—,—) sign 
structure for off-diagonal pairs. 

Predator-prey. The first entry in an off-diagonal pair 
is sampled from a half-normal distribution 1^(0, cr)| and 
the second entry from —|A/’(0, cr)| with probability C/2, 
or with the half-normal distributions reversed with prob¬ 
ability C/2, and both entries are zero otherwise. These 
community matrices have a (-1-,—) or (-,-!-) sign struc¬ 
ture for off-diagonal pairs. 


AY(t) = e’^‘AY(0) 
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C. Pseudospectra and transient instability 

In general, the eigenvalues of M satisfy the following 
definition: 

A(M) = {z e C : det(zl - M) = 0} (5) 

meaning that if z is an eigenvalue of M then by con¬ 
vention the norm of (zl — M)“^ is defined to be infinity 
(see Chapter /.I in [Ml)- But if ||(zl — M)“^|| is finite 
and very large, as is often the case with perturbed non¬ 
normal matrices, then the pseudospectrum of M must be 
considered. The ‘e-pseudospectrum’ has several equiva¬ 
lent definitions that describe the eigenvalues of a matrix 
whose entries have been subject to noise of magnitude 
e (in the sense of the matrix norm) [55]. We use the 
following definition: 

A,(M) = {z G C : ||(zI-M)-i|| > e-i} (6) 

If a matrix is normal then its e-pseudospectrum (hence¬ 
forth just ‘pseudospectrum’) consists of closed balls of 
radius e surrounding the original eigenvalues of M (see 
Theorem 2.2 in |29|1. As mentioned earlier, normal ma¬ 
trices cannot exhibit reactive dynamics: perturbations of 
the population vector for a stable system decay immedi¬ 
ately and with exponential profile as the system returns 
to its original equilibrium point. But with non-normal 
matrices, pseudospectra can be much larger and more in¬ 
tricate and reactive dynamics are possible: perturbations 
of the population vector for a stable system first increase 
in magnitude and reach a maximum amplitude before 
eventually decaying (Fig. 1). This behaviour motivates a 
description of local stability analysis for community ma¬ 
trices in terms of pseudospectra. 

Local asymptotic stability is determined in the same 
way for normal and non-normal matrices. The ‘spectral 
abscissa’ of M is defined as 

a(M) = sup Re(z) (7) 

zGA(M) 

where the supremum (sup) selects for the largest (real- 
part) of the rightmost eigenvalue in the set A(M). Sta¬ 
bility is guaranteed for Q!(M) <0. If M is normal, then 
||gMt|| _ gQ;(M)i dynamics are completely described 
by a(M) (see Eqn|^. Otherwise, the dynamics implied 
by M can be more complicated: 

ga(M)t < ||gMt|| < ^(V)e“(^)‘ (8) 

where the columns of matrix V are the eigenvectors of M, 
and «:(V) = ||V|| • ||V“^|| is known as the conditioning 
of V [55H5K] . The conditioning provides a bound from 
above—an upper bound—to the maximum amplitude of 
a perturbation of the population vector (it is worth not¬ 
ing that k(V) does not provide any information about 
the time at which the perturbation reaches its maximum 
amplitude). 

In complement to stability is reactivity, which de¬ 
scribes the behaviour of a system close to t = 0, at the 



t 

FIG. 1: Top: Pseudospectrum of a random community matrix 
with S = 50, C = 0.1, /r = 1 and a = 0.3, which is asymp¬ 
totically stable. Contours in the complex plane illustrate the 
effect on eigenvalues of the community matrix M for noise of 
magnitude e = 10’' m- The contour for e = 0.1 crosses the 
imaginary axis, implying that the pseudospectral abscissa is 
positive and so transient instability is observable. Bottom: 
Dynamics of ||e'^‘|| (arbitrary units of time, see Eqn[^. The 
dashed curve represents dynamics from eigenvalue analysis, 
whereas the solid curve represents dynamics predicted by pos¬ 
itive e-pseudospectral abscissa for e « 0.1. 


application of a perturbation. The ‘numerical abscissa’ 
of M is defined as 

w(M) = ^||e^*|| = sup Re(z) (9) 

at t=o 2 ga(h) 

where H = (n^TT]. The numerical abscissa is 

the maximum initial amplification rate following an in¬ 
finitesimally small perturbation to the population vector. 
Dynamics are non-reactive if a;(M) < 0 and may be reac¬ 
tive if a;(M) > 0. A stable system can be either reactive 
or non-reactive, but an unstable system is necessarily re¬ 
active. 

With non-normal matrices, perturbations to the en¬ 
tries of M can affect whether a system is stable and non¬ 
reactive or stable and reactive. In other words, perturba¬ 
tions to the entries of M can affect how a system responds 
to perturbations to the population vector. The effect of 
such perturbations to M is not covered by Eqn [^ How¬ 
ever, we can study the pseudospectrum of a community 
matrix to better understand system dynamics between 
the limits of reactivity and stability. In what follows, we 
use the theory of pseudospectra to relate uncertainty in, 
or small changes to, the entries of M to bounds on the 
amplification of perturbations of the population vector. 

The ‘e-pseudospectral abscissa’ of M is defined as 

a£(M) = sup Re(z) (10) 

2GA,(M) 
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which is the largest real-part eigenvalue of the pseu¬ 
dospectrum of M for a given amount of noise e. The 
e-pseudospectral abscissa provides a lower bound to the 
maximum amplification of a perturbation of the popula¬ 
tion vector (see Eqn 14.6 in [5^): 


sup^iiMl <sup||e^‘| 

£>0 e i>o 


( 11 ) 


and therefore the function 


/M(e) = 


a£(M) 


( 12 ) 


is useful for understanding transient dynamics [5S] . 
In the literature on pseudospectra, supj>g /M(e) = 
Aii(M) is known as the Kreiss constant |32l I34j. Eqns 
[TOl [TT] and \V2\ are useful because they relate perturbations 
to the matrix norm—small changes to the elements of the 
community matrix as described by the noise parameter 
e —to the effect of perturbations to the population vec¬ 
tor (compare Eqns and 10). For a given community 
matrix, as the size of a matrix perturbation is increased 
from zero there may be some critical value e* at which 
/M(e*) = 1. In the pseudospectrum, this is illustrated by 
the e*-contour crossing the imaginary axis (Fig. 1). At 
this point, perturbations to the equilibrium population 
vector begin to be amplified. 

For a stable and non-reactive system, perturbations to 
the population vector are not amplified and the system 
always returns to its original equilibrium point. For an 
unstable and necessarily reactive system, perturbations 
are amplified and the system may move to a new equilib¬ 
rium point. But for a stable and reactive system, pertur¬ 
bations are first amplified before the system eventually 
returns to its original equilibrium point—this is transient 
instability. Now that we can compute upper (Eqn|^ and 
lower bounds (Eqn 11) for amplifications, we are in a 
position to compare the transient dynamics of different 
types of ecological community described by non-normal 
community matrices. 


III. RESULTS 

We generated multiple sets of community matrices 
with C = 0.1, /r = 1 and various combinations of S and 
a for the five generative models. We first consider lower 
and upper bounds to the maximum amplitude of pertur¬ 
bations to the population vector for random community 
matrices, before turning our attention to the other types 
of interaction. The data required to reproduce the plots 
in this article are available at m- 


A. Lower bound for random community matrices 

We numerically evaluated the e-pseudospectral ab¬ 
scissa using the recently proposed subspace method [35] . 



FIG. 2: Regions of stability, transient instability and in¬ 
stability for a random community matrices with S = 100, 
C = 0.1 and /r = 1 as cr is varied. The y-axis is the lower 
bound of the maximum amplitude of perturbations to the 
population vector (Eqn [TT]). Transient instability is observ¬ 
able as the curve crosses one at an ~ 0.22 and instability 
is reached at a^ = ~ 0.31. At the threshold 

of instability, the lower bound of the maximum amplitude 
is lb((Jc) = 1.046 ± 0.006 (mean ± std). The shaded area 
represents the standard error over 100 realisations. 


Consider an ensemble of community matrices generated 
with random interaction type and S = 100 and a — 0.3, 
which is just below the threshold for instability {<7^ = 
-/= = « 0.31). We found that the average value 

of /M(e) (Eqn[l^ monotonically increases as a function 
of e and eventually saturates. At e* « 0.085 the curve 
crosses one, at which point perturbations are amplified 
and transient instability may be observable. The function 
/m (e) converges for all asymptotically stable community 
matrices considered here. 

In general, we identify regions of stability, tran¬ 
sient instability and instability by plotting supg>Q 
(Eqn 11 in practice, we plot fm{^) for large values of e) 
as cr is varied (Fig. 2). Similar regions can be identified as 
S is varied while a is held constant (results not shown). 
In the stable region, there is no perturbation to the com¬ 
munity matrix large enough (that can still be considered 
infinitesimally small) such that supg>Q > 1, and so 

perturbations are never amplified. At some critical point, 
CTti, there is a level of matrix noise e = e* above which 
perturbations to the population vector are amplified be¬ 
fore decaying. As a increases in the region of transient 
instability, e* decreases until it reaches zero at a^- At this 
point, system dynamics are guaranteed to be asymptoti¬ 
cally unstable and any infinitesimally small perturbation 
to the population vector is amplified (without necessarily 
returning to the original equilibrium point). In the un¬ 
stable region, /M(e) diverges and corresponding values 
for the lower bound should be treated with caution. 

The critical point for transient instability with S = 100 
is tJti ~ 0.22. This is very close to the value given 
by reactivity criteria based on the numerical abscissa: 
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FIG. 3: Distribution of upper bounds of the maximum am¬ 
plitude of perturbations to the popnlation vector (Eqn|^ for 
random community matrices generated with S = 100, C = 0.1 
and /r = 1 and seven values of cr (10,000 realisations). Distri¬ 
butions are fat-tailed and the slope of the tail does not change 
with (7. 


~ both approaches de¬ 

termine whether perturbations to the population vector 
are amplified based on eigenvalues related to M. As a 
point of difference, however, the pseudospectral approach 
allows for an additional treatment of uncertainty in, or 
small changes to, entries of the community matrix. For a 
given set of parameters, the numerical abscissa only indi¬ 
cates whether amplification is possible, whereas the pseu¬ 
dospectrum, through the e-pseudospectral abscissa, also 
indicates whether amplification is possible given small 
changes to the strengths of interactions among species in 
the community. 


B. Upper bound for random community matrices 

We plot the frequency distribution of k(V) (Eqn[^ for 
various combinations of S and a to investigate the upper 
bound to the maximum amplitude of perturbations of the 
population vector. In general, distributions are strongly 
peaked and fat-tailed (Fig. 3). This indicates that very 
large amplification is possible even for very small pertur¬ 
bations. The location of the peak changes very little as a 
increases, but shifts rightwards as S increases (results not 
shown). The slope of the tail does not change much as 
either S' or cr is varied. With S = 100 and a — ac = 0.31, 
the peak in the distribution of upper bound values is 
UBpeak(crc) ~ 95 and the maximum value in the tail is 
UBtaii(crc) 1000. When a power law is fit to the tail, 
f{x) oc x~°‘, the exponent is a « 2.9. 


TABLE I: Properties of community matrices with S = 100, 
C = 0.1, ^l=l. 


Type 

CTti 

CJ-c 

LB(crc) 

UBpeak(<^c^ 

UBtail c ; 

Q 

Mutualism 

0.11 

0.16 

1.02 

100 

~ 1000 

3 

Mixture 

0.17 

0.19 

1.02 

77 

~ 1000 

2.7 

Competition 

0.17 

0.20 

1.02 

100 

~ 1000 

3 

Random 

0.22 

0.31 

1.03 

95 

~ 1000 

2.9 

Predator-prey 

0.37 

0.87 

1.10 

60 

~ 500 

3.4 


C. Community matrices with different types of 
interaction 

The region of transient instability varies for different 
types of interaction, as do lower and upper bounds for 
amplification (Table 1). Transient instability becomes 
observable with smallest crti with mutualism, followed by 
mixture, competition, random and predator-prey. This 
order is the same as for the threshold for instability, Uc. 
However, the size of the region of transient instability, 
<7c — CTti, has a different order: predator-prey is largest, 
followed by random, mutualism, competition and mix¬ 
ture. The pattern is similar if S is varied while a is held 
constant (results not shown). As expected, these findings 
are consistent with earlier results based on the numerical 
abscissa and the correlation between off-diagonal entries 
in a community matrix m- 

Predator-prey community matrices are relatively sta¬ 
ble and exhibit the largest range of parameter values for 
transient instability. The lower bound to the maximum 
amplitude of perturbations of the population vector also 
reaches its largest value among the five types of inter¬ 
action for predator-prey community matrices. However, 
the peak in the distribution of upper bounds is at lower 
amplification and the slope of the tail is steeper (Table 1). 
This implies that perturbations are typically amplified 
less severely compared to the other types of interaction 
and the very largest possible amplitudes are not as large. 

Mutualism (+,+) and competition (—,—) have differ¬ 
ent critical points for transient instability and instability, 
but similar bounds to the maximum amplitude of pertur¬ 
bations of the population vector. Interestingly, the peak 
in the distribution of upper bounds is at lower amplihca- 
tion for community matrices with a mixture of these two 
interaction types. The largest upper bound, UBtaiKcTc), 
however, is similar to mutualism and competition, so the 
exponent, a, is shallower. 


IV. DISCUSSION 

Here we described transient instability for non-normal 
community matrices using local stability analysis and 
pseudospectra. We showed how the shift from stable and 
non-reactive dynamics to transient instability changes if 
perturbations are applied to the community matrix. We 
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also characterised how perturbations of the population 
vector are amplified during periods of transient instabil¬ 
ity for different types of interaction. We found an early, 
sharp and severe transition between stability and insta¬ 
bility with mutualism, mixture and competition, but a 
later, longer and less severe transition with predator-prey 
community matrices. 

In this study, we assumed a random topology of in¬ 
teractions between species. Although the correlation be¬ 
tween interaction strengths—and therefore the predom¬ 
inant type of interaction in a community matrix—may 
be more important than topology for stability nnumAt 
remains to be seen whether this is the case with transient 
instability. Nevertheless, it is likely that the particular 
trajectory of a perturbed system is sensitive to topology, 
and, of course, the direction of initial perturbation of the 
population vector. Understanding transient dynamics at 
this level of detail requires analysis of pseudoeigenvec¬ 
tors in addition to pseudoeigenvalues (see Chapter J.4 
in [H]). 

Local stability analysis is only one approach to under¬ 
standing the capacity for ecosystems to withstand exter¬ 
nal shocks [311 ED]- It will be informative to compare how 
the time evolution of the same shock to the same system 
is assessed under different approaches to measuring the 
‘stability’, ‘persistence’ or ‘resilience’ of ecosystems ini- 

stability, in principle, promises a degree of certainty 
that biodiversity will not be lost mm- Reactivity has 
been suggested as a possible early-warning signal for the 


onset of instability [TUIE5] . Transient instability not only 
fills the gap between these two concepts, but also high¬ 
lights new consequences of rapid environmental change. 
The longer the period of transient instability and the 
larger the amplification of perturbations of the popula¬ 
tion vector, the more susceptible an ecosystem is to mul¬ 
tiple perturbations. One perturbation may drive a stable 
system into a period of transient instability that eventu¬ 
ally dissipates; but two or three perturbations in quick 
succession may force the system to a new, unknown equi¬ 
librium point that may correspond to a loss of species 
and biodiversity. Pseudospectra can be used to investi¬ 
gate which ecosystems are at risk of instability, and what 
could be done to mitigate that risk. 
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